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Abstract 

The half filled Hubbard model is studied in the pair approximation of the Cluster 
Variation Method. The use of the SO (4) symmetry of the model makes possible to 
give a complete analytical characterization of the ground state, by means of explicit 
expressions for the double occupancy and the nearest neighbor correlation functions. 
The finite temperature analysis is reduced to the numerical solution of only two cou- 
pled transcendental equations. The behavior of local magnetic moment, specific heat 
and correlation functions is given for some typical cases in one and two dimensions. 
We obtain good qualitative agreement with exact and numerical results in one di- 
mension. The results for finite temperatures show a rapid evolution, with increasing 
temperature, from a strongly antiferromagnetic behavior to a disordered one; in the 
high temperature region a maximum (which has been related to a "gradual" metal- 
insulator transition) is found in the specific heat for very large values of the Coulomb 
repulsion. 
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1. Introduction 

The Hubbard model[l,2] is the simplest model of itinerant electrons which 
takes into account the interaction between electrons. It was originally proposed 
to describe the behavior of (i-electrons in transition metals, and it is expected 
to describe the metal-insulator (Mott) transition. In the recent years, the 
interest about this model has been greatly revived by the discovery of high-T c 
superconductors, since these materials are generally good Mott insulators and, 
in the superconducting phase, exhibit strong antiferromagnetic correlations, just 
like the half-filled Hubbard model at low temperatures. 

The model is defined by the following grand-canonical hamiltonian: 



where U, t > and a^, a ia and rii a are, respectively, annihilation, creation and 
number operators for electrons at site i with spin a G {+, — }. The first term 
represents the Coulomb repulsion between electrons at the same site (all other 
interactions are neglected); the second term is the chemical potential, and the 
third one is the kinetic term, which describes hopping of the electrons between 
sites, with the sum restricted to non-oriented nearest neighbor (n.n.) pairs. 

The Hubbard model has been studied by many different techniques (for 
reviews see[3,4] and references therein) but an exact solution is available only 
in one dimension [5, 6], while in two dimensions or more there are only a few exact 
results in very particular cases. For U/t = the model describes a system of 
non-interacting, moving electrons and is exactly solvable in any dimension. On 
the other side, for U/t = oo (atomic limit) and at half filling (i.e. (rii) = 
(rii + + rii-) = 1) the ground state is that of an antiferromagnetic insulator [7], 
with exactly one electron per site. At half filling two other very important 
rigorous results hold: 
(i) the chemical potential is given by ji = U/2 for any value of U/t and at any 

temperature [8] and 
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(ii) hamiltonian (1) has, for = U/2, an 50(4) symmetry [9]. 

In this paper we investigate the .D-dimensional Hubbard model at half fill- 
ing in the pair approximation of the Cluster Variation Method (CVM). The 
CVM has been originally introduced by Kikuchi[10] and its convergence in 
the thermodynamic limit has been demonstrated by Schlijper[ll]. Recently 
the method has been given a very elegant formulation by An[12], in terms of 
Mobius inversion. The simplest level of approximation in the CVM is the site 
approximation, which is equivalent to the ordinary mean-field theory; then we 
have the pair approximation, which can be shown to be equivalent to the Bethe 
approximation. 

The pair approximation of the CVM has already been applied to the 
Hubbard model in refs. [13-16]. Unfortunately, in these references only the 
U(l) ® U(l) Cartan subgroup of the SO (4) symmetry group later studied by 
Yang and Zhang was used, and the authors had to deal with large sets of cou- 
pled transcendental equations, which they could solve only for relatively high 
temperatures (kT/t > 1, with k Boltzmann's constant and T absolute temper- 
ature). Furthermore, in [14-16] equivalence is assumed between sites belonging 
to the two interpenetrating sublattices which form a bipartite lattice. 

In this paper we apply the pair approximation of the CVM to the Hubbard 
model on a bipartite lattice (that is, a lattice which is made of two interpenetrat- 
ing sublattices, say A and B, in such a way that a site belonging to sublattice 
A has all its nearest neighbors in the B sublattice and viceversa: examples 
are the 1-D chain, he, sq, sc and bec lattices) by taking into account the full 
SO {A) symmetry of the model. Furthermore, the equivalence between the two 
sublattices is not assumed, but it is derived from the thermodynamics of the 
model. The ground state is determined analytically, i.e. explicit expressions are 
derived for the double occupancy and the n.n. correlation functions at T = 0, 
for any number of dimensions and any value of the interaction U/t. For the 
finite temperature case, the problem is reduced to the numerical solution of 
two coupled transcendental equations. Such equations can be solved at any 
temperature and in the limit T — > the ground state solution is recovered. 

The validity of the approximation is first checked by comparing the behav- 
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ior of the zero temperature local magnetic moment vs. U/t in D = 1 with the 
exact solution for the infinite chain reported by Hirsch[17], and then by compar- 
ing the same quantity at finite temperature for typical values of U/t with the 
numerical results for finite chains obtained by Shiba and Pincus[18]. Once the 
validity of the method has been established, we report on the numerical results 
at finite temperature for D > 1: correlation functions, hopping expectations 
and specific heat are given in some typical cases, and the antiferromagnetic 
behavior of the system as well as the inhibition of certain hopping processes at 
low temperatures are discussed. Finally, it is noticed that a high temperature 
maximum appears in the specific heat for very large values of the interaction, in 
agreement with previous studies where this maximum was related to a gradual 
metal-insulator transition. 

The paper is organized as follows: in Sec. 2 we construct the trial free 
energy according to the pair approximation of the CVM, taking into account 
the symmetry of the hamiltonian. In Sec. 3 the ground state is obtained and 
discussed, and the zero temperature local magnetic moment is compared with 
the exact result for the infinite chain. In Sec. 4 the analysis is extended to 
finite temperature and the behavior of various physical quantities is given and 
discussed and, finally, in Sec. 5 some conclusions are drawn. 

2. Free energy 

Following An's formulation[12], the CVM trial free energy for a bipartite 
lattice in the pair approximation can be written as 

/ = | + k B T^-^- {^{pAlrvpA) + Tr(p B lnp B )] + ^Tr(p p ]np p ) } , (2) 

where E is the internal energy, N is the number of lattice sites, and pa, Pb and 
p p are the reduced density matrices (to be determined by minimizing /) for a 
site belonging to sublattice A, a site belonging to sublattice B and a pair of 
nearest neighbors, respectively. 

Before taking the variation of / with respect to the reduced density ma- 
trices, let us determine which constraints for such density matrices can be de- 
rived from the symmetry group of the hamiltonian. As shown in [9], hamilto- 
nian (1), with p = U/2 because of the half filling condition, commutes with a 
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SO (A) = — group, where one SU(2) (referred to as the magnetic 

Z-<2 

one) is generated by 

Jz = 2 J+ =^24 + ai_, J_ = ^a\_a i+ (5) 

i i i 

the other SU(2) (called the pairing, or superconductive one), which is relevant 
only at half-filling, by 

2 : J>i + +n*--l), K+ = J2e^al + at, K_ = e** i a i -a i+ (6) 

i i i 

(the phase factor e 1 ^ is +1 for sites in sublattice A and —1 for sites in sub- 
lattice B) and TL-i interchanges the two SU(2) symmetries. The presence or 
absence of this symmetry in the quantum state of the system characterizes the 
different phases: a disordered phase will be invariant under the whole 50(4) 
symmetry group, while a phase with magnetic and/or superconductive order 
will be invariant under a reduced symmetry group, with both SU (2) (or one) 
spontaneously broken down. Since in this paper we are concerned with the half 
filled case, and it is conjectured that at half filling the Hubbard model does not 
undergo any phase transition, we devote our attention to the disordered phase, 
thus assuming the whole SO (4) symmetry (in view of the good agreement with 
exact results in one dimension, this assumption should be correct also at zero 
temperature). The possibility of a phase transition, associated with a sponta- 
neous breaking of the magnetic SU (2) symmetry group, will be examined in a 
forthcoming paper[19] for the extended Hubbard model at general filling. 

In order to impose the commutation relations between the reduced density 
matrices and the SO (4) generators defined above, we introduce in the site and 
pair reduced Fock spaces the customary basis of eigenstates of the number 
operators. In such a basis, requiring that the reduced density matrices commute 
with the Cartan operators J z and K z , p& and ps turn out to be diagonal, 
while for p p one obtains the same block structure as in [16], with only 36 non- 
zero elements. By imposing furthermore the commutation with J + and K + 
(or, equivalently, with their hermitian conjugates J_ and K_), one finds that 

p 7 (7 = A, B) has two distinct eigenvalues, say <i 7 and <i 7 , each with 

multiplicity 2, whereas p p is a block diagonal matrix with: 
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i) two eigenvalues Ai and A2, each with multiplicity 2; 



ii) four degenerate 2x2 blocks, which give rise to two eigenvalues A3 and A4 
with multiplicity 4; 

iii) a 4 x 4 block with eigenvalues Ai, A 2 , A 5 and A 6 . 

Summarizing, we have 6 different eigenvalues, Ai and A2 with multiplicity 
mi = m 2 = 3, A3 and A4 with multiplicity m 3 = = 4, and A5 and A6 with 
multiplicity = rriQ = 1. Of these, only 5 are independent, because of the 
normalization condition Tr(p p ) = 1. 

Recalling that the expectation value of an operator X is given by (X) = 
Tr(pX), one can compute, with some simple algebra, the expectation values of 
all the site and n.n. pair operators. The non-zero expectation values turn out 
to be (i and j nearest neighbors, i E A and j G B) 



(n ia ) = (n jo ) = - 
(m+rii-) = d A 
(n j+ nj-) = d B 

c p = (n i+ n j+ ) = (ni-rijJ) = Ai + A 2 + A 3 + A 4 



c a = (m+rij-) = (n j+ m-) = Ai - A 2 H 

(n i+ ni-n j+ ) = (m+m-nj-) = ^ (d A + c p + c a - ^ 

{n j+ nj-n i+ ) = (n j+ 7ij-7ii-) = \ (d B + c p + c a -\ 



1 - d A - d B 



(7) 





q = (n i+ ni-n j+ nj-) = X 1 



for the diagonal operators (notice that in this scheme the half filling condition 
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is derived from the symmetry and not imposed), and 



P = {a\ + a\_aj-a j+ ) = ^ - c p - c a 
p' = {a\_a] + a i+ aj-) = c a - c p 

TT = (aL^ni-aUj-a) = {a\ a a ja {l - n l _ fJ )(l - n 3 - a )) 



1 



= \\J(^-^) 2 -\(d A -d B y 



(8) 



y = ( a L a jM ~ ni- a )nj- a } = (a^a^m-ail - n J _ (T )) 
= \ \A A s - A 6 ) 2 - [2(d A + d B )-l- 3(Ai - A 2 )] 2 

(together with the obvious hermitian conjugates) for the non-diagonal opera- 
tors. 

The free energy per site, as a function of d A , d B and Aj, i = 1, . . . 6, is then 



/ = | (d A + d B -l)- 2zt^ (A 3 - A 4 ) 2 - \(d A - d B y 
- ztyj (A 5 - A 6 ) 2 - [2(d A + d B )-l- 3(Ai - A 2 )] 2 



+ fcT(l-z) £ 

7 =A,B 
6 

+ feT| ^(miAi hi A;) 

i=l 



<i 7 In <i 7 + ( - — c2 7 J In f - — d- 



(9) 



3. The ground state 

At T = 0, the free energy per site is but the internal energy, and is given 



by 



/ = ^(d A + d B - 1) - 2ztJ (A 3 - A 4 ) 2 - - rfs) 2 

2 y 4 (io) 

- ^ (A 5 - A 6 ) 2 - [2(d A + d B ) - 1 - 3(A X - A 2 )] 2 . 

The ground state can thus be obtained by minimizing / with respect to the d 7 
and the Aj . Since we are looking for an absolute minimum and our variables are 
subject to constraints (the eigenvalues of the density matrices, as well as the 
arguments of the square roots in (10) must be non-negative, and p p must be 



properly normalized), we should search our minimum possibly at the domain 
boundary of the constrained variables, and not only in the interior. 
Indeed, the minimum is found for 

^ K 1 - vm) • " = S< 

A5 = 1 and Aj = 0, % 7^ 5. 

The ground state is thus described by 

1 / U 

d>A = d,B = d = — 1 — 



4 V VW 2 + 16 



1,1/. W 
= - -rf= - 1 + 



2 4 V VU 2 + 16/ (12) 

ri = ^2^(1 - 2d) = 



VK 2 + 16 
9 = c p = r = 0. 

The configuration of a pair of nearest neighbors in the ground state can be 
derived as the eigenvector of p p corresponding to the eigenvalue A5 = 1. One 
obtains, up to a normalization constant, 

2 



* = 



. ia\ + a\_ + a],a]_) 

VK 2 + IQ + 3+3 

If. U " 



|0). 



VW 2 + 16, 

It is worth noticing that such a configuration is the superposition of an 
antiferromagnetic singlet pair (equivalent to that used by Anderson [20] in the 
construction of his RVB state) with a fraction of doubly occupied sites. The 
concentration of doubly occupied sites as well as the kinetic energy (the expec- 
tation value of the hopping term) vanish for large U and have their maximum 
for U = 0. The n.n. correlations are strictly antiferromagnetic (i.e., c p = 0), as 
expected, and no phase transition is found in any number of dimensions. Some 
hopping processes, e.g. the hopping of an electron from a doubly occupied to a 
singly occupied site, or from a singly occupied to an empty site, turn out to be 
inhibited in the ground state. 

In order to check the validity of our approximation, we compare the local 
magnetic moment at T = vs. U/t for the 1-D chain (z = 2) with the exact 
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result reported in [18]. The local magnetic moment S is proportional to the 
expectation value of the square of the magnetization: 



5=|((n i+ -n i _) 2 > (14) 
and is directly related to the double occupancy, since 

S = 7(1 - 2d) = I ( 1 + , f ^ . (15) 
4 V 8 V VU 2 + 16/ 

(15) is exact both in the non-interacting case (W = 0, £ = 3/8) and in the 
atomic limit (U = 00, S = 3/4). Fig. 1 shows the comparison between our 
results for z = 2 (solid line) and the exact solution for the 1-D chain (circles). 
The agreement is within 10% for all values of U/t. 



4. Finite temperature 

Since we have found cIa = ds = d in the ground state, and the entropy 
contribution favors this latter condition, one can expect this symmetry relation 
to hold even at finite temperature. Furthermore, a breaking of such symme- 
try at finite temperature would yield a reentrant phase with staggered double 
occupancy, and there is no indication of such phases in the Hubbard model. 
Indeed, we have checked numerically that the minima of / always appear for 
dA = ds = d. We shall therefore assume from now on the latter relation. 

In this way we obtain a free energy which is a function of six indepen- 
dent variables only: d and five of the A^. Instead of minimizing / directly, we 
introduce the following new set of indipendent variables: 

S = 4d-1 

r n , n+ i = m n (X n - A n+i ), n = l,3, 5 (16) 

Rn,n+1 = "^n(A n + A n+ i), U = 1, 3, 

and we define = 1 — R12 — -R34. After rewriting /, as given by (2), in terms 
of the above variables, the minimum-/ requirement gives (assuming, with no 
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loss of generality, A 3 > A 4 ) 









dr 12 2n 4 R 12 

or 34 2 4 R 34 - r 34 



ri2 



df 
dr 56 

df 
dR 12 

df 



zt ^ kT- In ~^ 56 ~*~ r56 
2ri 4 .R 56 - r 56 



(17) 



In fl?2 



12 



36 



-In 



56 



56 



34 



j n -^34 r 34 _ j n -^56 r 56 

64 " " ' 4 



Upon defining a; = T\ 2 — 5 and after some algebra three of the above equations 
can be solved for r 34 , R12, -R34 and 5 (or n.2) leaving us with the following two 
coupled transcendental equations for x and r^: 



x = tanh 



Z (3 
2(z-lf\n 2 



e.Rsinh (3— 



x 



r 56 = 2i2sinh I /?— 

' Ti 



(18) 



where /5 = (kT/t) 1 , T\ = -\Jt%q — x 2 and 



Tl 



6 cosh (3— + 8 cosh (3 + 2 cosh /3 



,T 5 6 
Tl 



(19) 



Once (18) has been solved, the remaining variables are given by the following 
relations: 



7*12 = — 6-Rsinh ( (3 — 



X 



R12 = G.Rcosh 



x 



(20) 



r 3 4 = 8R sinh (3 
R 34 = Stfcosh/?. 

As a check for the whole procedure, we compare in Fig. 2 our results for the 
local magnetic moment for z = 2 and some typical values of U/t with the exact 



(numerical) results obtained by Shiba and Pincus[18] for a six sites chain with 
periodic boundary conditions. We find good qualitative agreement, and again 
differences are contained within 10%. It can be observed that the solution for 
low temperatures converges to the value predicted by the ground state analysis. 
The results for the chain are compared with those for the square lattice (z = 4) 
in Fig. 3: our analysis shows, as expected from numerical simulation[17], that 
the local moment increases with increasing U/t and with decreasing dimension- 
ality. 

In Fig. 4 we report the correlation functions c p (lower curves) and c a (upper 
curves), in Fig. 5 the hopping contributions tq (lower curves) and t\ (upper 
curves), in Fig. 6 the "double hoppings" p (lower curves) and p' (upper curves) 
and in Fig. 7 the specific heat, for U/t = 8 (solid lines) and U/t = 4 (dashed 
lines) for a square lattice. 

Of course there is no evidence of a true phase transition (the specific heat 
exhibits a maximum but not a sharp peak) , but we can clearly distinguish a low 
temperature behavior (kT/t < 0.5) from a high temperature one (kT/t > 1). 
The low temperature region is characterized by strong antiferromagnetic corre- 
lations and by a relatively large kinetic energy associated to the moving elec- 
trons, due almost entirely to double hoppings and to hopping processes from 
doubly occupied to empty sites and viceversa, while the remaining processes 
are strongly inhibited because of the ground state configuration. The high tem- 
perature region, besides, looks like a true disordered phase, with almost equally 
distributed correlations (c p ~ c a ~ 1/4) and low kinetic energy. Furthermore 
in this region, as already noticed in [16], for very large values of the interaction 
U/t a spread maximum appears in the specific heat, which was related by Ho 
and Barry to a "gradual" metal-insulator transition. 

5. Conclusions 

We have investigated the half-filled Hubbard model in the pair approxima- 
tion of the Cluster Variation Method, making use of the full SO (4) symmetry 
of the model. We have given an analytical description of the ground state, by 
means of the double occupancy and of the n.n. correlation functions and, for 
finite temperature, we have derived a pair of coupled transcendental equations. 
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Numerical solution shows two different behaviors, connected by a smooth but 
rapid change in the values of the parameters. The low temperature behavior is 
strongly antiferromagnetic and exhibits the inhibition of certain hopping pro- 
cesses, while a large kinetic energy is associated with the others. In the high 
temperature region we find a quite disordered behavior, with a spread maxi- 
mum, which was related to a metal-insulator transition, for very large values 
of the interaction. Good agreement is found with exact and numerical results 
in one dimension. 
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Figure Captions 

Fig. 1: Local magnetic moment at T = 0. Our result (solid line) and result 
from [17] (circles). 

Fig. 2: Local magnetic moment at finite temperature for U/t = 8 (upper 

curves) and U/t = 4 (lower curves). Solid lines are our results, 

circles are from [18]. 
Fig. 3: Local magnetic moment at finite temperature for U/t = 8 (upper 

curves) and U/t = 4 (lower curves). Dashed lines are for the linear 

chain and solid lines for the square lattice. 
Fig. 4: Correlation functions c p (lower curves) and c a (upper curves) on the 

square lattice for U/t = 8 (solid lines) and U/t = 4 (dashed lines). 
Fig. 5: Hopping expectation values To (lower curves) and T\ (upper curves) 

on the square lattice for U/t = 8 (solid lines) and U/t = 4. 
Fig. 6: Double hoppings p (lower curves) and p' (upper curves) on the square 

lattice for U/t = 8 (solid lines) and U/t = 4 (dashed lines). 
Fig. 7: Specific heat on the square lattice for U/t = 8 (solid lines) and 

U/t = 4. 
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